% Estimate dynamic factor as described in Ritschl, Sarferaz, and Uebele (2015)
% Constant version

clear all; clc;

addpath('../funcs');
addpath('../scrpts');
addpath('../data');
%% Load and prepare data
name53_1867_1995_new

subset = 'all'; % choose subset
large=0; % choose dataset (for large=1 set subset = 'data98')

data53 = load('data/US_53_1867_2006_version_8_2009.txt');
data98 = load('data/US_98_1867_1939_version_10_2008.txt');

if large==1
    data=data98;
else
    data = data53;
end
% Apply HP filter
data_hp = dataPrep(data,6.25,15,2,3);

% Define subsets
setsubset_new;
%subset = 'all';

if size(data_hp,2)>53
    yt = data_hp';
    
    filename = char('matFiles/constant/results_US_98_1867_1929_constant');
else
    yt = data_hp(:,eval(subset))';
    
    filename = char('matFiles/constant/results_53_1867_2006_constant_');
    filename=char([filename,subset]);
end


%% Gibbs sampling preliminaries
chain           = 100000;   % number of draws for gibbs sampler
burn            =  chain*0.8; % burn-in phase
jumpsave        = ceil((chain-burn)/1000);    % number of draws that are skipped
chain_save      = (chain-burn)/jumpsave;
gx              = 1;

%% Model specifications
K = 1; % number of factors
q=8; % lag length of factors
p=1; % lag length of idiosyncratic components
m = K*q;

sv = 0;
stat = 0; % stationarity check (stat=1) and without (stat=0)

%% Number of observations
[N,T] = size(yt);
n = K*(p+1);
%T = T+1;
%% Priors

% Prior for the VAR parameters
tau_f = 0.2; % uncertainty around prior mean of zero of factor coeff.
tau_u = 1; % uncertainty around prior mean of zero of idio. coeff.

% Prior for (V)AR coefficients

% Zero restrictions
VPhi_prior = ones((K)*q,K)*1e-12;

% Prior for those parameters without zero restrictions

% Prior for Phi
for ix=1:K
    for qx=1:q
        VPhi_prior(K*(qx-1)+ix,ix) = tau_f/qx;
    end
end

% Prior for Theta
%for ix=1:N
for px=1:p
    VTheta_prior((px-1)+ix,ix) = tau_u/px;
end
%end


Phi_prior = zeros((K)^2*q,1);
VPhi_prior = diag(vec(VPhi_prior));

Theta_prior = zeros(p,1);
VTheta_prior = diag(vec(VTheta_prior));

Lam_prior = zeros(N*K,1);
VLam_prior = eye(N*K,N*K)*100;

Tu_0= 6;% prior for shape of VCV
Qu_0=eye(N)*0.001;%*(Tu_0-p-1);% prior for scale of VCV


%% Initial conditions
f0 = zeros(m,1);
Pf0= eye(m);

%% Starting Values
ft_pc = extract(yt')';
[maxcorr,indx] = max(corr(ft_pc',yt'));

ft    = ft_pc;

%lamt   = randn(N,K,T);

Phi = VAROLS(ft',q,1);
Theta = eye(N)*0.1;
Omega_chi = eye(N)*0.1;


%% Pre-allocate Space

results.ft          = zeros(K,T,chain_save);
results.lam        = zeros(N,K,chain_save);
results.Phi         = zeros(K,K*q,chain_save);
results.Theta       = zeros(N,N*p,chain_save);
results.Omega_chi   = zeros(N,N,chain_save);
results.nut         = zeros(N,T-1,chain_save);

results.stato_f     = zeros(1,chain_save);
results.wsign       = zeros(1,chain_save);

lamvec = zeros(N*K,1);
tic
for w=1:chain
    
    if mod(w,100)==0
        [ssff errorm]   = sprintf('Sampler at %d out of %d draws, subset: %s',w,chain,subset);
        disp(ssff);
    end
    %% Quasi-difference data
    yt_star = quasiDiff(yt,T,N,p,Theta);
    
    %% Draw factor loadings
    for ix =1:N
        ft_star = ft(2:end) - Theta(ix,ix)*ft(1:end-1);
        
        if ix==1
            lamvec(ix,:) =   multiReg(yt_star(ix,2:end)',ft_star',Omega_chi(ix,ix),VLam_prior(ix,ix),Lam_prior(ix),2);
        else
            lamvec(ix,:) =   multiReg(yt_star(ix,2:end)',ft_star',Omega_chi(ix,ix),VLam_prior(ix,ix),Lam_prior(ix),2);
        end
    end
    
    lam = reshape(lamvec,[N,K]);
    
    %% Draw factors
    ft = filterSmootherCK(yt_star,T,N,m,K,q,f0,Pf0,@state_space_f,ones(1,T),Omega_chi,Phi,Theta,q,p,repmat(lam,[1 1 T]));
    
    
    %% Draw (V)AR coefficient of factors
    [Phi, stato_f, nu_star] = alphaVAR(ft',zeros(size(ft')),K,q,VPhi_prior,Phi_prior,stat,sv);
    
    %% Draw (V)AR coefficient of idionsyncratic components
    ut = zeros(N,T-1);
    for tx=1:T-1
        ut(:,tx) = yt(:,tx)-lam*ft(:,tx);
    end
    
    for nx=1:N
        Theta(nx,nx) = alphaVAR(ut(nx,:)',Omega_chi(nx,nx),1,p,VTheta_prior,Theta_prior,stat,2);
    end
    
    
    %% Draw variance of idionsyncratic equation
    Omega_chi = alphaVariance(ut,Theta,eye(N),p,Qu_0,Tu_0);
    
    %% Save draws
    
    if w>burn
        if mod(w,jumpsave)==0
            results.ft(:,:,gx) = ft;
            results.lam(:,:,gx) = lam;
            results.Phi(:,:,gx) = Phi;
            results.Theta(:,:,gx) = Theta;
            results.Omega_chi(:,:,gx) = Omega_chi;
            
            gx                      = gx + 1;
        end
    end
    
end

toc
save(filename ,'results', 'yt','T','q','N','K','Phi_prior','VPhi_prior',...
    'Theta_prior', 'VTheta_prior', 'Qu_0', 'Tu_0',...
    'f0', 'Pf0','sv', 'ft_pc', 'chain', 'chain_save');

